function [hfile, densm] = srf_ToSphere(hfile, niter, mforce, sforce)
% SRF::ToSphere  - perform to-sphere morphing
%
% FORMAT:       [srf, densm] = srf.ToSphere([niter, mforce, sforce]])
%
% Input fields:
%
%       niter       number of iterations, default [12000, 8000, 5000, 3000]
%       mforce      morphing force, default [0.5, 0.25, 0.1625, 0.125]
%       sforce      to-sphere force, default [0.0003, 0.0025, 0.01, 0.02]
%                   arrays must match in size if given
%
% Output fields:
%
%       srf         spherical surface
%       densm       density SMP
%
% Note: this method simply passes to SRF::Morph. It automatically set
%       area and distortion correction (to 3 - eps)

% Version:  v0.7b
% Build:    7093002
% Date:     Sep-30 2007, 2:49 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% global settings from config
global bvqxfile_config;

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
ngiven = false;
if nargin < 2 || ...
    isempty(niter) || ...
    numel(niter) > 8 || ...
   ~isa(niter, 'double') || ...
    any(isnan(niter(:)) | niter(:) < 0 | niter(:) > 1e5)
    try
        niter = bvqxfile_config.settings.Morphing.ToSphere.NrOfIterations;
    catch
        niter = [12000, 8000, 5000, 3000];
    end
else
    ngiven = true;
    niter = floor(niter(:)');
end
if nargin < 3 || ...
   ~isa(mforce, 'double') || ...
    numel(mforce) ~= numel(niter) || ...
    any(isnan(mforce(:)) | mforce(:) <= 0 | mforce(:) >= 1)
    if ngiven
        mforce = 0.5 * (niter:-1:1) ./ niter;
    else
        try
            mforce = bvqxfile_config.settings.Morphing.ToSphere.Force;
        catch
            mforce = [0.5, 0.25, 0.1625, 0.125];
        end
    end
end
if nargin < 4 || ...
   ~isa(sforce, 'double') || ...
    numel(sforce) ~= numel(niter) || ...
    any(isnan(sforce(:)) | sforce(:) <= 0 | sforce(:) >= 1)
    if ngiven
        sforce = 0.02 * (2 ^ 0:(niter - 1)) ./ (2 ^ (niter - 1));
    else
        try
            sforce = bvqxfile_config.settings.Morphing.ToSphere.SphereForce;
        catch
            sforce = [0.0003, 0.0025, 0.01, 0.02];
        end
    end
end
opts = struct;
opts.areac = 1;
opts.distc = 3 - eps;
title = 'Morphing to sphere';

% get filename
sc = bvqxfile_getscont(hfile.L);
srffile = sc.F;
if numel(srffile) > 4 && ...
    strcmpi(srffile(end-3:end), '.srf')
    tfilename = [srffile(1:end-4) '_ToSphere.srf'];
else
    tfilename = '';
end
% iterate over steps
for nc = 1:numel(niter)
    opts.sphere = sforce(nc);
    opts.title = sprintf('%s, forces: %f/%f', title, mforce(nc), sforce(nc));
    if nc < numel(niter) || nargout < 2
        srf_Morph(hfile, niter(nc), mforce(nc), 'smooth', opts);
    else
    [hfile, densm] = srf_Morph(hfile, niter(end), mforce(end), 'smooth', opts);
    end
    if ~isempty(tfilename)
        try
            aft_SaveAs(hfile, tfilename);
        catch
            % do nothing
        end
    end
end

% center sphere and force to radius
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
sc.F = srffile;
c = bc.VertexCoordinate;
mnc = min(c);
mxc = max(c);
c = c - repmat((mnc + mxc) / 2, [size(c, 1), 1]);
cr = sqrt(sum(c .* c, 2));
mcr = mean(cr);
c = c .* repmat(mcr ./ cr, [1, 3]);
c = c + repmat(bc.MeshCenter, [size(c, 1), 1]);
sc.C.VertexCoordinate = c;
bvqxfile_setscont(hfile.L, sc);
srf_RecalcNormals(hfile);

% remove tempfile
if ~isempty(tfilename) && ...
    exist(tfilename, 'file') == 2
    try
        delete(tfilename);
    catch
        % do nothing
    end
end
